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ABSTRACT 



Two atmospheric prediction models, based upon the meteorological 
primitive equations, were programmed and tested on the CDC 6500 computer, 

A frictionless, barotropic model was integrated, using ten-minute time steps, 
at 500 MBS for periods up to four days. Then, an n -level - , baroclinic model 
for an inviscid, adiabatic atmosphere was integrated for brief test periods at 
five uniformly-spaced zeta (sigma) surfaces over a smooth earth. Both models 
were initialized with actual data fields provided by FNWC Monterey. 

Acceptability was based on measurements of energy and vorticity para- 
meters , as well as on qualitative assessments of output fields. The 
barotropic 500 MB height forecasts were found to be comparable to the FNWC 
(vorticity model) barotropic forecasts. Mean-square -vorticity and kinetic 
energy were suitably conserved for integrations up to four days. An energy 
accumulation in the smaller range of scale was increasingly evident beyond 
two days into the forecast period. No smoothing was performed to control 
this accumulation, A limited number of test integrations was made with 
the five-level model. Computational instabilities were observed for 
forecasts beyond twelve hours. Interim results are presented. 
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1.0 INTRODUCTION 



Disadvantages arising from more generalized forms of the vorticity, 
divergence, and thermodynamic equations have overshadowed the advant- 
age of relatively large integration time steps in filtered prediction models. 
At the same time, increased knowledge in finite -difference techniques, 
in initialization procedures, and in boundary condition formulation, have 
been realized in a decade marked by unbelievable increases in computer 
capacity and speeds. These trends and accomplishments have accelerated 
the shift to the so-called "primitive equations" as a basis for atmospheric 
prediction models, both operational and research. 

1.1 Finite Differencing 

Whereas, in filtered models, an integration time step of one hour is 
possible without risk of linear computational instability, the corresponding 
time step for primitive -equation integrations is approximately ten minutes. 
One must fulfill the von Neumann stability criterion in order to avoid this 
type of instability in linear differential equations . A second type of 
instability may arise in spite of careful selection of time step and grid 
mesh length. This is the non-linear computational instability, or 
"aliasing", first described by Phillips (1959), Although Phillips treated 
non-linear differential equations, Miyakoda (1962) showed that non-linear 
instability can also occur in special types of linear equations. 

Aliasing is described by Lilly (1965) for regions in which some arbi- 
trary continuous function is being represented by a finite number of 
gridded values. It is shown that there exists a maximum and minimum 
wavelength which can be resolved by any grid array, and that waves 
outside of this range will be misrepresented in terms of the permitted 
harmonics . 
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Phillips controlled this phenomenon by smoothing. It may be reduced 
significantly, however, by use of conservation forms of the difference 
analogs, especially those which preserve the average wave number. The 
latter prevent the cascade of energy toward smaller ranges of scale. The 
well-known methods of Arakawa (1966) represent those which conserve 
mean-square vorticity, total kinetic energy, and average wave number. 
Phase errors are not eliminated by conservation forms, but their 
magnitudes are tolerable. 

Another potential difficulty is associated with the phrasing of time 
derivatives, which have taken on a new significance with the advances 
made in space -differencing methods. The work of Kurihara (1965), Lilly 
(1965), and Young (1968), is illustrative of efforts to study, and possibly, 
to reduce the "slow 1 ' instability associated with the various time -differenc- 
ing schemes. The "leapfrog" scheme of Richtmyer (1963) has been used 
widely in spite of its shortcomings. These may be listed as follows: 

a. kinetic energy is not conserved, 

b. during extended integrations, non-linear instability may develop 
in the absence of some form of smoothing. 

c. the odd- and even-solution sets tend to proceed independently of 
each other because of the alternating sign of the computational mode. 

(some form of solution averaging may be helpful to eliminate the latter 
difficulty). In spite of these shortcomings, the leapfrog scheme is used 
because of the implied economies and its general accuracy. 

The finite-differencing scheme ultimately selected for a given model 
is one which best compromises its associated: (1) diffusive character- 
istics; (2) existence and properties of any computational mode; (3) 
stability properties; and (4), the computational burden. 
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1.2 Boundary Conditions 

The boundary conditions for any model must be tailored both to the 
specific equation set and to the objective of the integrations (short range 
forecast, or long-term general circulation experiment). The primitive- 
equation free -surface barotropic forecast model of Shuman and Vanderman 
(1966) is indicative of the effects of approximate versus correct boundary 
conditions on the forecast fields. Shuman showed that the RMS velocity 
divergence was tolerably conserved for periods up to about thirteen days 
when using approximate boundary conditions, but that the model "blew 
up" rapidly therafter. By using the correct boundary conditions (which 
were consistent with the equations) , the model not only conserved RMS 
velocity divergence from the outset, but also prevented further deterioration 
of the forecast which used the approximate conditions for the first thirteen 
days . 

The important point is that incorrect boundary conditions (as well as 
discretization errors) may excite gravity -inertial waves which, in the 
absence of some dissipative mechanism, may contaminate the forecast. 

See Nitta and Hovermale (1967) for an interesting discussion of these 
matters. Matsuno (1966) reported on false reflections of the computation- 
al mode from an outflow boundary. This inward propagation of errors was 
especially prominent at the shorter wavelengths. The rate of reflection 
decreased with increasing wavelength and with increasing "order" of 
continuity at the lateral boundaries. Thus, to overcome this effect, an 
increasing number of interior points is required to approximate the 
boundary value of an advected parameter (if a suitable assumption or 
approximation cannot be found for the particular problem). 
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1.3 Initialization Procedures 

Initialization procedures are probably the least well-known of the 
requirements for model development- The whole problem of weather 
analysis and forecasting is concerned with the mutual adjustment between 
the mass and motion fields under the influence of rotation and gravity. 

One very interesting result arose from the general circulation experi- 
ments of Smagorinsky (1965). Patterns of precipitation appeared roughly 
the same regardless of the initial condition of vertical motion, provided 
that other initial conditions (such as temperature and the non-divergent 
wind component) are held fixed. Miyakoda and Moyer (1968) suggest 
that for periods of more than five days (for the integrations) , the vertical 
motion initialization is somewhat arbitrary because of the deterministic 
nature of the equations. 

Initialization falls into two classes: (1) extrapolative, and (2) 
iterative. In either case, the idea is to reduce the noise and to represent 
only the slowly changing meteorological state. The former type consists 
of a marching integration of the equations to arrive at this "balanced" 
condition. The latter type stays at the same time level, but iterates 
the computations until the equilibrium state is reached. Clearly, the 
choice depends upon the model purpose. The iterative method would 
have an advantage for users requiring accurate short-term forecasts (of 
the order of a few hours, say), whereas the extrapolative method would 
be entirely satisfactory for forecast periods beyond the time required 
to "settle down" (ten hours or longer, say). 

Thus far, primitive -equation models have been used with an assort- 
ment of initial fields - some real and some analytical constructions. 
Linearized research models may be initialized according to some pre- 
selected analytical distribution with a known solution. 
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In contrast, Shuman and Hovermale (1967) chose to use as initial data 
the objective analyses of geopotential height and temperature at the ten 
mandatory pressure levels (from 1000 to 100 MBS). Each of the input 
analyses was expanded from the JNWP octagonal grid into a 53 x 57 
rectangular array by an extrapolative scheme which consisted of: (1) 
finding the mean value of the analyzed quantity on the octagon boundary, 
and copying this value into the non-overlapping portion (border region) 
of the 53 x 57 array; and (2) "wedding" the two regions by means of a 
"ring" smoother, while holding the boundary of the rectangular array 
fixed. By this process, one recognizes that the equatorial regions 
would be virtually devoid of small scale components. 

The removal, or reduction, of these small-scale features is clearly 
justified for several reasons, not the least of which is the general 
inappropriate ness of the so-called balance equation in the tropics. The 
difficulties, moreover, of false reflections as reported by Matsuno (1966) 
would be minimized in such "flat" fields. By flattening a field in this 
way, one clearly "increases the order of continuity" at the expense of 
accurate specification of meteorological detail in the tropics in order 
to preclude possible computational difficulties. 

The finite -difference schemes, boundary conditions, and initializa- 
tion procedures for the models studied in this paper will be discussed 
in detail in the appropriate sections of this thesis. 

1.4 Coordinate Systems 

Historically, the (x,y,p,t) coordinate system was shown to be 
superior to the rectangular system in many respects, but it suffered 
from the difficulties arising from the lower boundary condition. 
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Over smooth terrain, one generally assumes that the vertical velocity 
vanishes identically at the lower boundary; whereas, over rough terrain, 
one assumes that 



co = \V • ^7 -rr 



(or some equivalent approximation) where IT is the terrain pressure. 

In Phillip's (x,y,<T,t) sigma system, the particular problem just 
mentioned is overcome, but the sigma system is not entirely free of 
practical difficulties. Kurihara (1968) reported, for example, on problems 
arising from (large) terrain-pressure gradients. 

In this research project, the barotropic model was written in 
pressure coordinates, but the baroclinic model was based upon a 
generalized "zeta" system which was proposed by R.T. Williams 
(personal communication). In this modified sigma system, the vertical 
coordinate zeta is related to Phillip's sigma by an arbitrary functional 
relationship 



where, it is clear, one may achieve the desired "linearity" in the 
vertical parameter while allowing for arbitrary placement of the 
computational surfaces . 
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2.0 THE BAROTROPIC EXPERIMENT 



2.1 Introduction 

The matter of testing a simple barotropic model prior to engaging in 
more complex simulations is a desirable approach. The work of Houghton 
et al (1966) and Shuman and Vanderman (1966) is cited in this regard. 

Houghton conducted experiments with various finite -difference 
formulations for the Coriolis term and boundary conditions in order to 
observe sensitivities in long-term solutions. Solution stability was 
achieved without a large amount of smoothing, but some was necessary 
in order to control truncation errors. Charney (1965) had suggested an 
upper limit (several weeks) for long-term integrations because of 
inaccuracies in initial conditions. Houghton allowed that stability and 
truncation effects may be of equal importance to the validity of long- 
term solutions as to the growth of initial error. 

In this project, the author has speculated that the small errors, or 
inconsistencies, in the initial conditions, as well as the effects of 
truncation and discretication, will not be of major concern for short- 
term integrations. There does appear to be some basis for more concern 
about the former in one- and two-day forecasts. 

Shuman and Vanderman (1966) also tested a primitive -equation free- 
surface barotropic model in a tropical channel experiment. They 
concluded that approximate boundary conditions (which were inconsistent 
with the difference analogs) were entirely satisfactory for integrations 
up to about twelve days. The philosophy behind Shuman’s choice was 
that the conditions provided for reflection of pure gravity waves at the 
walls. The inertial effects (due to the presence of "f"), it was reasoned, 
could be ignored for a limited time. 
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The finite -difference analog referred to by Shuman as the "semi- 
momentum" form apparently contains a large amount of smoothing 
because of the generous use of the superbars. This may account no, 
only for the success of the long-term integrations without explicit 
smoothing, bu, also for the degree of concern about "unsmoothing" 
before use of the balance equation. 



2.2 Purpose of the Barotropic Experiment 

The first phase of thesis research consisted of the derivation, 
programming (CDC 6500), and evaluation of a barotropic free-surface 
model. This was test run at the 500 MB level using, as initial data, 
the actual objectively -analyzed hemispheric data fields provided by 
FNWC Monterey. 

The author, after noting the many achievements of investigators 
cited earlier, felt that the obstacles remaining for reasonably-accurate 
short-term integrations of real data fields were more practical than 
theoretical. If would instructive, that is, to see if approximations 
and assumptions found suitable in analytical test integrations were 

equally suitable for integrations with real data fields . 

The purpose of the barotropic experiment was to note the suitability 
of certain boundary conditions, to determine appropriate initialization 
procedures, and to test the modified-Arakawa finite -difference analogs 
(as proposed by R.T. Williams) for a model which might be expected to 
reasonably simulate atmospheric flow patterns at 500 MBS for periods 
up to four days . 
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2.3 The System of Prognostic Equations 



The two-dimensional differential equations of momentum and 
continuity for a free-surface barotropic, inviscid atmosphere may be 
approximated on a polar stereographic projection as follows: 

a. momentum equation in east-west direction: 



a L,vU __ 
£-fc 






^ < d 
- rm 






5 I )2 — \ v_ o> / ] 

rm ) ao V ,wt ) J 






d 3 



(2.3.1) 



b. momentum equation in north-south direction: 



h 5 + (W>a + ur(u ^ + IT ^ 

1 <) <^j a* ag 

tdx\ /vy\ } dg\ am ) ) 




(2.3.2) 



c. continuity equation: 

~ - (YA \ _^L al/\) 4- / 

Ux a 3 v ) 



(2.3.3) 



The reader is referred to Appendix A for the development of this set 
of equations in flux form. 

The computational boundary was taken to be mid-way between the 
outer ring of grid points and the first interior ring of points . This 
permitted the practical application of boundary conditions required for 
conserving the square of any parameter being advected, while at the 
same time, permitting non-zero normal wind components at the first 
interior ring. 
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Geostrophic flow was assumed for the normal wind component along 
the first interior ring of points because of the computational convenience. 
By so doing, the problem of computing the correct values of "h" along 



the outermost ring is avoided. 

The advective terms led to the following boundary conditions . At 
the computational boundaries parallel to the y-axis , 




(2.3.4) 



while, at the boundaries parallel to the x-axis. 




(2.3.5) 



For convenience, the tangential wind component at the first interior 
ring was copied into the adjacent grid points on the outermost ring. In 
a similar manner, the values of "h" were also shifted "out" for 
convenience . 

It will be seen that since the minimum value of the Coriolis 
parameter permitted was that value corresponding to "f" at 15° latitude, 
the map factor advection terms (which are small anyway) vanish south 
of that latitude. The local change of "uh" and "vh" depended, therefore, 
on the horizontal flux divergence terms at the first interior ring (since 
geostrophic flow was assumed). Single (space) differences were used 
as appropriate in the continuity equation along the first interior ring. 
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Initialization of wind components consisted of solving the so-called 



linear balance equation of the form 



V <§ = 7 * (+7+) 



(2.3.6) 



using an FNWC 500 MB height field. The "first guess" streamfunction 
values were taken to be 



where h was the initial height value at each point. 

The integration time step was ten minutes. Subsequent to an 
initial forward (single) step, a conventional centered difference was 
used. 

2.4 Finite Difference Forms 

A modified-Arakawa finite -difference scheme was employed. In 
this scheme, the momentum and continuity equation analogs were 
expressed in the forms shown below: 

a. momentum equation in east-west direction: 




(2.3.7) 



o 






(2.4.1) 
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b. momentum equation in north-south direction: 







where, for any advected parameter, 

2 » 



t 



we have defined the operator 





where "n" denotes the time level. 






( 2 , 4 , 4 ) 
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In an attempt to reduce some of the small scale irregularities of the 
initial height analyses near the computational boundary, a five-point 
smoother routine was written to produce the following effect. First, a 
mean height value (say) was calculated for the outermost grid ring. 

This value was set at all of the outermost ring points. The smoothing 
coefficient was latitude dependent in the sense that the smoothing 
increased from zero at 20° latitude to a degree sufficient to remove 
all perturbations at the outermost ring. This procedure was deemed to 
be consistent with, and justified by, the linear balance assumption 
in the tropics . 

2.5 Research Procedures and Results 

The barotropic model was tested on actual FNWC 500 MB height 
analyses for the Northern Hemisphere, rather than on some idealized 
distribution for which some analytical solution was known. This less 
rigorous test was consistent with the purpose of the study - to see if 
some of the approximations found suitable for analytical tests were 
equally suitable for use with actual data fields (and their much larger 
frequency ranges). 

The judgement as to solution stability or acceptability was based 
on certain diagnostic parameters, such as the square -vorticity 
parameter (SVP) or mean kinetic energy (KE) , as well as on the meteor- 
ological appearance of the output fields. 

The first test consisted of computing all parameters and terms as 
they appeared in the equations (without space averaging). Then, the 
model was rerun using space averaging in a manner similar to that 
used by Shuman (1966). It was determined that the SVP grew at an 
intolerably large rate. 
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When no space averaging was performed, the SVP growth rate was an 
order of magnitude smaller. In spite of the lesser rate, however, energy 
did tend to accumulate in the smaller range of scale, particularly in 
tropical regions . 

A second test consisted of measuring the SVP for a run in which the 
integration cycle was restarted every 72 time steps, but with no space 
averaging. This was an attempt to measure the effects of non-linear 
interactions between the odd- and even-solution sets. The result was 
that the SVP grew at approximately the same rate for both runs, but the 
SVP took a small step upward in value immediately following the 
restart. 

Figure 1 depicts the aforementioned SVP growth rates. The upper 
curve (A) shows the rate for the case when space averaging was used, 
but no restarting was performed. The lower curve (B) shows the SVP 
rate for both cases in which no space averaging was performed (no 
attempt was made to indicate the small upward jogs after each restart, 
since they were negligible on this scale). The larger rate (A) was 
judged to be the result of some rather serious imbalances caused by 
not "unsmoothing" initial fields in the manner reported by Shuman (1966, 
1967). The much smaller rate (B) was most likely caused by small 
imbalances arising from the use of the linear balance equation, and by 
the gradual accumulation (which could be either a real or computationally- 
induced cascade) of energy in the small scale features. No smoothing 
was performed in case B. 
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Figure 1. Square-Vorticity Parameter (SVP) as a function of time step. Curve A pertains to 
case with space averaging of terms. Curve B applies to cases in which no space averaging 
was employed; one in which the cycle was restarted after 36 steps; and one in which the 
cycle was never restarted. No significant SVP difference occurred in the latter cases. 




Figure 2. Difference Field (in meters) between the odd- and even- 
solution sets after 144 ten-minute time steps in barotropic experi- 
ment. This area shown represents a lOd x lOd subset of the total 
,?mispheric field. It is idealized to portray the typical scale 
and magnitude of such a difference field. In general, the range of 
sc^li of features is from 2d to 4d, with magnitudes generally less 
than five meters. 
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A third test was made to measure the "departure" of the two 
solution sets. Figure 2 illustrates the (idealized) scale and magnitudes 
of the difference field after 144 time steps. Each side of the figure is 
ten mesh lengths long (about 3810 kilometers) , and is representative of 
the entire hemispheric difference pattern. In general, the magnitudes 
of the centers are five meters or less after twenty-four hours of forecast 
period. This may be compared to the magnitude of the 500 MB diurnal 
height oscillation, which is (roughly) an order of magnitude larger than 
this difference field. 

The mean kinetic energy was computed for each time step. No 
attempt was made to plot these distributions because the variation was 
negligible in the cases in which no space averaging was performed. 

One test was attempted in which the latest odd- and even-step 
wind component and height fields were averaged (as appropriate) after 
72 time steps, and then restarted. The calculations, however, blew 
up immediately after restarting the cycle. Thus, it would appear as 
though the "average" wind components and "average" heights 
introduce a considerable trauma in the system, in general. It should 
be noted that, at all other times, the restarting procedure involved 
using the latest fields rather than averaged fields. 

Charts A through E depict the initial 500 MB analysis and daily 
forecast fields out to four days, respectively. The model tends to 
overdevelop both highs and lows , but it is felt that the phenomenon 
could be easily controlled. Note also that the amplitudes of small 
scale features increase perceptibly (with time), particularly in low 
latitudes. No attempt was made to control this accumulation of energy. 



29 



The barotropic primitive -equation 24-hour 500 MB forecast fields 
were compared with the barotropic filtered model output fields provided 
by FNWC . The height error fields revealed that the PE model output 
compared favorably with the operational FNWC model in the case 
tested. Charts F and G illustrate the height error patterns for the 
PE and filtered models, respectively, for the 24-hour period ending 
0000Z, 25 April 1968. It was beyond the intent of this research to make 
more than a perfunctory comparison of these two models. 

In summary, it can be said that the initialization through use of the 
linear balance equation was satisfactory when small-scale features 
were removed from the boundary regions by a latitudinally-varying 
filter. The boundary conditions and finite -difference forms used in 
the experiment were both economical and stable for integrations up to 
four days. Although the smaller-scale features tended to increase in 
amplitude, it was felt that the phenomenon could be easily controlled. 
The results would justify further research on a barotropic model 
initialized with real hemispheric data fields , 
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Chart A. Initial 500 MB Height Analysis for 0000Z, 24 April 1968. 
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CHART B 




Chart B. 



P.E. Barotropic 500 MB 24-Hour Height Forecast from 
0000Z, 24 April 1968. 
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Chart C . 



P.E. Barotropic 500 MB 48-Hour Height Forecast from 
0000Z, 24 April 1968. 
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Chart D. P.E. Barotropic 500 MB 72-Hour Height Forecast from 
0000Z, 24 April 1968. 
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Chart E. P.E. Barotropic 500 MB 96-Hour Height Forecast from 
0000Z, 24 April 1968. 
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Chart F . 



P.E. Barotropic 24-Hour 500 MB Height Error Field 
(in meters) for period ending 0000Z, 25 April 1968. 
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Chart G. FNWC (Vorticity Model) Barotropic 24-Hour 500 MB 
Error Field (in meters) for period ending 0000Z, 

25 April 1968. 
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3,0 THE BAROCLINIC EXPERIMENT 



3.1 Introduction 

A generalized Fortran program was written for integrating an n-level 
baroclinic primitive -equation set for a frictionless, adiabatic atmos- 
phere. Because of its generality, much of the possible computational 
economy was sacrificed. As a result, the frequent use of auxiliary 
memory (disc file) caused the expenditure of about twice as much 
peripheral processor time as central processor time per time step . 

To check out and test run the model, the input-output section 
was written for five uniformly-spaced zeta (sigma) surfaces only. The 
vertical coordinate (zeta) may be given by the transformation 



where <T is defined as the ratio of pressure to that at the lower 
boundary (IT); that is, 




(3.1.1) 




(3.1.2) 



Differentiation with respect to sigma leads to 




(3.1.3) 



which, in turn, suggested the definition of the 
map factor" 



so-called "vertical 



This leads to the derivative transformation 



ah 







as 



( 3 . 1 . 5 ) 



where A is an arbitrary dependent variable. In the special case used 
for evaluating the model, the transformation was linear and the vertical 
map factor, s, was equal to unity. Moreover, since zeta increased 
upward, the vertical velocity was defined positive upward according 
to the relation 



It is well-known that Phillip’s sigma system (and therefore, the 
zeta system) has both advantages and disadvantages. Its main 
advantages are: 

a. the absence of one of the chief difficulties in pressure 
coordinates - that of needing to define the horizontal pressure gradient 
at sea level. For models that include terrain, the generation of 
ficticious pressure gradients is avoided because the terrain pressure 

is observed rather than calculated through some reduction scheme. 

b. the function representing vertical motion vanishes identically 
at all material surfaces. 

Two of the disadvantages of the sigma system are: 

a. that observational data are generally taken at pressure 
surfaces. Thus, sigma surface parameters must be generated by some 
interpolation/extrapolation scheme . 

b. that computational difficulties may arise because of large 
space variations of terrain pressure in certain regions (see Kurihara, 
1968 ). 



UT 3 - cr 



( 3 . 1 . 6 ) 
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3.2 Purpose of the Baroclinic Experiment 

Realizing that time would not permit the inclusion of diabatic heating 
and diffusion terms in the equations, the goal was made to program 
an n-level model for an adiabatic, inviscid atmosphere on a smooth 
earth; and to test the model at five uniformly-spaced zeta surfaces (in 
the time remaining). Provisions were made in the Fortran program to 
incorporate the remaining terms later without need for any major 
reprogramming effort. 

As in the barotropic experiment, it remained to continue the 
evaluation of boundary conditions and initialization procedures. The 
modified-Arakawa conservation-type difference scheme had to be tested 
in the multi-level version. In addition, the author was aware of the 
buildup of energy in small-scale features in the barotropic model, and 
desired to observe the manifestations of frictionless flow in the more 
complex model, 

3.3 The System of Prognostic Equations 

The equation set for an inviscid, adiabatic atmosphere may be 
approximated on a polar stereographic map projection as follows: 
a, east-west momentum equation: 

Bitu _ _ iY M a /TTu. 



+wU-< ( £^ y) 
l 2^ 



w(rr £$ + RT 



c>TT 

gK 



(3.3.1) 
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b. north -south momentum equation: 



^EiTs: - (r»M 4 -TT5 

( dy\ rm ) JJ ^ 






c?3 



c. continuity equation: 

d. hydrostatic equation: 

a§ = rt 

s<r 

e. thermodynamic equation: 

(Mu # (S &) f - 



. SI f-UV-j-C*! rWu ^IT + iT 

ay a$ 



>1 



where all symbols have been defined previously. 



(3.3.2) 



(3.3.3) 



(3.3.4) 



(3.3.5) 



Figure 3 contains a diagram of computational levels with notations 
as to which variables are predicted or computed at those levels. 
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Level 

Computed Variables (Notational 



Zeta 

1-0 r 


(at each level) 


Subscript) 


Sigma 

- n.n 


0.9 — . 


u, v,T, z 


5 




n , a _ _ 


w 


4 


n ? 


0.7 - 


u, v,T, z 


4 




n r _ 


w 


3 


n 4 


0.5 - 


u, v,T, z 


3 




0.4 


w 


2 


- n 5 


0.3 


u, v,T, z 


2 


0.7 


0.2 


w 


1 


n q 


0.1 


U, V,T, z 


1 


0.9 


0.0 


Boundary Pressure (pi) 







Figure 3. Diagram Zeta (Sigma) Levels. Notations are made to 
indicate the subscripts used for the level, as well as to show the 
variables predicted or computed at each level. Note that "z" has 
been used interchanaeably with phi (geopotential), and that "pi" 
ha ; been used for the 1 iwer boundary pressure. 
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3.4 Finite -Difference Forms 

The difference equations corresponding to (3.3) are based on the 
Arakawa technique. It will be noted that the set (3.3) are similar in 
form to those used earlier by Smagorinsky (1965). 

The difference scheme avoids non-linear instability by requiring 
that the advective terms conserve the square of the advected property 
(when summed over the entire grid), assuming continuous time 
derivatives. The total energy is conserved because of requirements 
placed upon the vertical differencing; specifically, upon the special 
form of the hydrostatic equation. Assuming that the zeta surfaces are 
uniformly-spaced, the lowest phi distribution is given by 

& = $ 4- ( -X-') (3.4.1a) 

' ~SFC ^ K S<3* M 



and for subsequent levels, one uses 




It should be noted that a more complicated form of the hydrostatic 
equation is required if the surfaces are not uniformly-spaced. The 
reader is referred to Haltiner (3) for a complete discussion of the 
conservation properties of this difference scheme. 

The east-west momentum equation is given by 



n'H-i hi - 1 

(ITU) r + 

fr ■** 




( irx-uv )»1 

2a 1 - J 





where (*) denotes the point (i,j,k). 
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Similarly, the north-south momentum equation is written 



-fe^i! [<$ - $ V (IT; .. -It. . ) l 

The lower boundary pressure is computed from 



( 3 . 4 . 3 ) 



M 



/Vv4t M-l 

\1 ‘j ^ l J 



( 3 . 4 . 4 ) 



where, for the entire column, one evaluates the local change as 
<V\ 5 ^ % ^v\ 



(52) a -S 5<x -o< +,4 -|S 

at .j A ^ L *hjk 7 jH 

where, for notational convenience, one defines the quantities 



( 3 . 4 . 5 ) 



and 



ck - 


U1T 

Av\ 


/3- 


l/j? 

/W\ 



The thermodynamic equation takes on the form 

ttfT) = (ifT)^ + a it 



kP 



"M- + ^ 'Vj J 



/V\ 



( 3 . 4 . 6 ) 
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Note that for any advected parameter an operator has been defined 

• 2 . 

[ik+. + )- 



+[L ft + ,+fyXfy .+ty - ( ) J . , 



A 'id > t j 



(3.4.7) 



Finally, vertical motions for non-material surfaces are obtained 
from a form of the continuity equation 

ur„ = ur _ $ /5t\ ^ 

* SJ' r IW'q 

j 

~ (3.4.8) 

. 

3.5 The Computational Sequence 

For the initialization portion of the program, the sequence of steps 
is as follows: 

a. compute map factor and Coriolis parameter. 

b. read input fields from magnetic tape: 

(1) sea level pressure. 

(2) temperatures at 850, 700, 500, 400, 300, 200, and 

100 MBS. 



c. convert temperatures . 

d. interpolate p-surface temperatures to build zeta-surface 
temperatures, and smooth fields in tropics. 
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e. form T-pi fields, 

f . integrate hydrostatically to obtain phi fields , 

g. solve linear balance equation at each level, 

h. compute u- and v-components for each level, 

i. form u-pi and v-pi fields. 

For the steps involved in each time step after initialization is 
completed, the following are listed: 

a. compute local change in lower boundary pressure. 

b. compute vertical motion fields. 

c. compute next lower boundary pressure, 

d. compute forcing functions for momentum and thermodynamic 
equations for this level. 

e. integrate to obtain next u-pi, v-pi, and T-pi fields for 
this level, 

f. form next u, v, and T fields for this level. 

g. repeat steps d. through f. for remaining levels. 

h. compute next phi fields, 

i. (optional) compute diagnostic parameters. 

It should be noted that the integration cycle was restarted with a 
single forward time step every six hours in the baroclinic experiment, 

3.6 Boundary Conditions 

The boundary placement is identical to that used in the barotropic 
experiment. As in the barotropic experiment, the advective terms led 
to the following boundary conditions . 
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At the computational boundaries parallel to the y-axis , 



* 




(3.6.1a) 



and 







ax 



c? 



while, at the boundaries parallel to the x-axis, 







(3.6.1b) 



(3.6.2a) 



(3.6.2b) 



The practical imposition of boundary conditions was varied slightly 
from one test to another. For the most part, however, geostrophic 
flow was assumed at the first interior ring in order to avoid the 
necessity of estimating the values of "phi" and "pi" along the outer- 
most ring. One-sided space derivitives were used in the continuity 
and thermodynamic equations (as appropriate) for computations along 
the first interior ring. 

At the earth's surface and at the top of the atmosphere, the vertical 



velocity, w, vanished identically. 



3 . 7 Initialization Procedures 



The input fields for the baroclinic experiment were actual data 
fields as analyzed on the FNWC 63 x 63 hemispheric grid. To remove 
most of the disturbance components from the boundary regions , a 
latitudinally -dependent smoother of the form 

S(4>)=- < 3 - 7 -» 



was employed to compute the coefficient. The coefficient vanished 
north of 30° latitude. 

Initialization of the geopotential and wind fields was accomplished 
by use of the linear balance equation for a sigma(zeta) surface, (It 
should be noted that the terms arising from the presence of the 
sphericity terms in the momentum equation were not retained.) The 
balance equation took on the form 

= jr(j§ + RT WW 



+ -LVT.n7£*tt 
T 









(3.7,2) 



as a result of the divergence of the two-termed pressure gradient force 
in sigma coordinates. 

Appendix C contains a derivation of the divergence equation in 
sigma coordinates. 
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3.8 Research Procedures and Results 



The generality of the Fortran program led to a very inefficient use 
of the CDC 6500 computer at the Fleet Numerical Weather Central. For 
the five-level version tested, approximately 180 disc accesses per time 
step were required (for storage of the history variables, and temporary 
storage of intermediate results). Because of the extremely heavy 
burden of high-priority operational programs being run on the FNWC 
computer in the period, October 1968 to December 1968, this phase 
of research was limited. In addition to "debugging" the 2000-instruction 
program, a few integrations were run out to twelve hours in order to 
make a few assessments of the output fields and measurements of the 
conservation parameters . 

Interim results of the baroclinic experiment indicate the generation 
of static instabilities in the uppermost layers, and a subsequent rapid 
deterioration in the fields. No attempt was made to constrain the 
lapse rates of predicted temperatures. 

The output fields at the sixth hour into the forecast did not show 
any indication of computational deterioration. 

In almost every instance, the computational blow-up was 
accompanied by a catastrophic potential-to-kinetic energy conversion 
in the uppermost layers . 

A continuing series of test runs are being made by the author in 
an effort to ascertain the source of the computational difficulties. 
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4.0 CONCLUSIONS AND RECOMMENDATIONS 



4.1 The Barotropic Experiment 

The barotropic experiment was intended to test the suitability of 
boundary conditions, initialization procedures, and the phrasing of 
the finite -difference analogs. In addition, some information as to 
the degree of smoothing of the initial data fields in the tropical regions 
was sought. 

Initialization through use of the linear balance equation was 
satisfactory when most of the disturbance component was removed 
at the boundary region through use of a latitudinally-varying smoother 
routine . 

The boundary conditions and finite -difference forms were found 
to be both economical and stable for integrations up to four days . (No 
runs were made beyond four days.) 

Although energy tended to accumulate in the smaller ranges of 
scale, it was felt that this phenomenon could be controlled. 

The results of this investigation would justify further research on 
a barotropic primitive -equation model that was initialized with real 
hemispheric data fields. The conservation properties of the difference 
scheme place the model at an advantage over the operational FNWC 
model for intermediate length integrations (three to seven days). 

4.2 The Baroclinic Experiment 

Since most of the test runs made indicated that instabilities 
developed in the first eight to thirteen hours of the forecast period, 
more work is required to ascertain the nature and cause (s) of the 
computational instabilities. 
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Preliminary results would seem to indicate that the vertical velocity 
fields were not satisfactory, generally . This was felt to be particularly 
true in the upper atmosphere, where sudden transformations of potential 
to kinetic energy were linked to unrealistic vertical transports, and the 
ensuing static instabilities. 

The generalized Fortran program should be useful as a vehicle for 
future graduate research work at the Naval Postgraduate School. 
Provisions have been made for arbitrarily varying the number and 
location of zeta levels. In addition, heating and friction terms may 
be added without need for extensive reprogramming. 

An extended core system (ECS) was added to the CDC 6500 computer 
in late December 1968, and will provided an additional 500,000 words 
of fast-access auxiliary memory. With little additional programming, 
this should greatly reduce the "field handling" time, and permit longer 
test runs . 
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APPENDIX A 



Flux Forms of Barotropic Equations in Pressure Coordinates 



The flux forms of the barotropic primitive -equation set in pressure 
coordinates may be "transformed" from the standard forms for a polar 
stereographic projection. These forms to be transformed are given by: 




(A .la) 



(A. lb) 



(A. 2) 



The procedure consists of showing that local change terms may be 
written as: 



5a = _L 5 _ u ^ , 

^ l 5-t 'j-f- j 



— 

5fc 



k l 



IT 




(A. 3a) 



(A. 3b) 



Similarly, the advective terms lead to the forms, 

W ^£4 L T - '111 j \ / k g) I Ivmt \ 

[\ \_ \ Sx\iw / ^ { rw\ \ 

4 - Jri \ 4“ ^( u lt 2vy)\1 

rtY)^ l <^4r ^ >) '] y 



(A. 4a) 
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and 



^ <3g l, 



(jwjf \ j^_ \ 

3 x V / 3 1_] \ ‘V\-\ ) 



4- 



iT ) 



Ml 



ib + k( 



u ^ + u-^vn 



2* 






^ "1 

>5 



(A. 4b) 



These terms may be substituted, as appropriate, in (A.l) to obtain: 









u + ir ^JV»\ 1 1 

a 5 'J 



-rmM 

( ^ v /vv> ) ^ 



(A. 5a) 



•^7p' l = -U)+n 4- (wits + cr(u _l (7 <3wa\1 
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where the continuity equation is given by. 



(A. 5b) 



r - (YA> 



> -f-(wU') + A-(LrU) } 



(A. 6) 



It should be noted that (A. 6) differs slightly from the usual flux form 
in which the map factor appears in both flux quantities. Thus, the 
initial divergence will not vanish identically (after use of the balance 
equation) «. 
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APPENDIX B 



Energy Computations 



In pressure coordinates the total potential energy (internal plus 
geopotential) may be written (for a smooth earth) as, 

rh 




CpT dp dA 



(B.l) 



AO 

where "dA" represents an area increment. Now, if one eliminates 
from the hydrostatic equations in pressure and zeta coordinates, 
respectively, one derives the relationship, 




(B.2) 



where "s" is the vertical mapping factor, and “TT is the lower boundary 
pressure (sea level). Substituting (B.2) in the first expression, the 
total potential energy becomes , 

r _ Si f (' U d? dA (B.3) 

f a ) ) S 

A o 

In a similar manner, the appropriate form for kinetic energy becomes, 

Ek= H drdA (B,4) 

A o 
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The corresponding finite -difference forms for (B , 3) and (B.4) were 




where "N" is the number of mesh areas and " L" is 
layers. In addition, "D" is the grid mesh length, 
were as initially specified. 



the number of 
and other symbols 
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APPENDIX C 



Divergence Equation Forms in the Sigma System 

One can take the divergence of the vector equation of frictionless, 
adiabatic motion on a sigma surface and show that the resulting 
divergence equation assumes the form, 

+ ‘S 1 - -tST - a o(u,w) + vo - . - u if 

ci-fc f ^<T <5X 




+ 







if the "sphericity" terms are omitted from the equation of motion. 

It has been shown (see Haltiner (3), for example) that the terms of 
significance to less -generalized vorticity models specialize to 



a 

Vf = 




(C.2) 



usually referred to as the linear balance equation, under the assumption 
that the non-divergent component of the wind is large compared to the 
divergent component. 

In sigma coordinates, however, terms arising from the additional 
pressure force term in the momentum equations are rather large, and 
must be included in the balance equation in sigma coordinates. The 
resulting form is given as 

<?<f --f-Vf V ktIv-Shtt + -Lyr-Vi.ifl' 

L- T 
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If one applies the (l^‘ Vx ) operator to the same equation of 
motion, the vertical component of the relative vorticity on a sigma 



surface is found to be 

£ - i- + KT -j 



S \ ^ - 



±i¥T.VJU 



(c . 4) 



which preserves the mathematical form of the linear balance equation. 
It should be noted that (C . 3) takes on a more complicated form when 
the term 



( try - U.V) 
Za 



(u,v) 



is added to the momentum equations . 



* as appropriate 
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beyond twelve hours. Interim results are presented. 
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